<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01//EN"
   "http://www.w3.org/TR/html4/strict.dtd">

<html>
<head>
  <title></title>
  <meta http-equiv="content-type" content="text/html; charset=latin1">
  <style type="text/css">
td.linenos { background-color: #f0f0f0; padding-right: 10px; }
span.lineno { background-color: #f0f0f0; padding: 0 5px 0 5px; }
pre { line-height: 125%; }
body .hll { background-color: #ffffcc }
body  { background: #f8f8f8; }
body .c { color: #008800; font-style: italic } /* Comment */
body .err { border: 1px solid #FF0000 } /* Error */
body .k { color: #AA22FF; font-weight: bold } /* Keyword */
body .o { color: #666666 } /* Operator */
body .cm { color: #008800; font-style: italic } /* Comment.Multiline */
body .cp { color: #008800 } /* Comment.Preproc */
body .c1 { color: #008800; font-style: italic } /* Comment.Single */
body .cs { color: #008800; font-weight: bold } /* Comment.Special */
body .gd { color: #A00000 } /* Generic.Deleted */
body .ge { font-style: italic } /* Generic.Emph */
body .gr { color: #FF0000 } /* Generic.Error */
body .gh { color: #000080; font-weight: bold } /* Generic.Heading */
body .gi { color: #00A000 } /* Generic.Inserted */
body .go { color: #888888 } /* Generic.Output */
body .gp { color: #000080; font-weight: bold } /* Generic.Prompt */
body .gs { font-weight: bold } /* Generic.Strong */
body .gu { color: #800080; font-weight: bold } /* Generic.Subheading */
body .gt { color: #0044DD } /* Generic.Traceback */
body .kc { color: #AA22FF; font-weight: bold } /* Keyword.Constant */
body .kd { color: #AA22FF; font-weight: bold } /* Keyword.Declaration */
body .kn { color: #AA22FF; font-weight: bold } /* Keyword.Namespace */
body .kp { color: #AA22FF } /* Keyword.Pseudo */
body .kr { color: #AA22FF; font-weight: bold } /* Keyword.Reserved */
body .kt { color: #00BB00; font-weight: bold } /* Keyword.Type */
body .m { color: #666666 } /* Literal.Number */
body .s { color: #BB4444 } /* Literal.String */
body .na { color: #BB4444 } /* Name.Attribute */
body .nb { color: #AA22FF } /* Name.Builtin */
body .nc { color: #0000FF } /* Name.Class */
body .no { color: #880000 } /* Name.Constant */
body .nd { color: #AA22FF } /* Name.Decorator */
body .ni { color: #999999; font-weight: bold } /* Name.Entity */
body .ne { color: #D2413A; font-weight: bold } /* Name.Exception */
body .nf { color: #00A000 } /* Name.Function */
body .nl { color: #A0A000 } /* Name.Label */
body .nn { color: #0000FF; font-weight: bold } /* Name.Namespace */
body .nt { color: #008000; font-weight: bold } /* Name.Tag */
body .nv { color: #B8860B } /* Name.Variable */
body .ow { color: #AA22FF; font-weight: bold } /* Operator.Word */
body .w { color: #bbbbbb } /* Text.Whitespace */
body .mf { color: #666666 } /* Literal.Number.Float */
body .mh { color: #666666 } /* Literal.Number.Hex */
body .mi { color: #666666 } /* Literal.Number.Integer */
body .mo { color: #666666 } /* Literal.Number.Oct */
body .sb { color: #BB4444 } /* Literal.String.Backtick */
body .sc { color: #BB4444 } /* Literal.String.Char */
body .sd { color: #BB4444; font-style: italic } /* Literal.String.Doc */
body .s2 { color: #BB4444 } /* Literal.String.Double */
body .se { color: #BB6622; font-weight: bold } /* Literal.String.Escape */
body .sh { color: #BB4444 } /* Literal.String.Heredoc */
body .si { color: #BB6688; font-weight: bold } /* Literal.String.Interpol */
body .sx { color: #008000 } /* Literal.String.Other */
body .sr { color: #BB6688 } /* Literal.String.Regex */
body .s1 { color: #BB4444 } /* Literal.String.Single */
body .ss { color: #B8860B } /* Literal.String.Symbol */
body .bp { color: #AA22FF } /* Name.Builtin.Pseudo */
body .vc { color: #B8860B } /* Name.Variable.Class */
body .vg { color: #B8860B } /* Name.Variable.Global */
body .vi { color: #B8860B } /* Name.Variable.Instance */
body .il { color: #666666 } /* Literal.Number.Integer.Long */

  </style>
</head>
<body>
<h2></h2>

<div class="highlight"><pre><span class="kn">from</span> <span class="nn">libatomism</span> <span class="kn">import</span> <span class="o">*</span>
<span class="kn">from</span> <span class="nn">math</span> <span class="kn">import</span> <span class="o">*</span>
<span class="kn">import</span> <span class="nn">array</span> 
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="kn">as</span> <span class="nn">plt</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="sd">This script shows how to model a molecule, the hydrogen molecule is taken as example. </span>
<span class="sd">The following feature are described:</span>
<span class="sd">    - managment of units</span>
<span class="sd">    - z-matrix generalized coordinates</span>
<span class="sd">    - equation of motion: kinetic operator and potential energy surface</span>

<span class="sd">We recomand you to use the interpreter ipython for the numerous sugar features it offers.</span>

<span class="sd">Some tips:</span>
<span class="sd">    - all objects can be printed by &#39;print myobjetc&#39;</span>
<span class="sd">    - to obtain the documentation of a class, input &#39;myClass?&#39; and press enter (ipython only)</span>
<span class="sd">    - to obtain the list of method in an object, input &#39;myobject.&#39; and press tab (ipython only)</span>
<span class="sd">    - to obtain the doxumentation on a method in an object, input &#39;myobject.myMethod?&#39; and press enter (ipython only)</span>
<span class="sd">&quot;&quot;&quot;</span>


<span class="sd">&quot;&quot;&quot; Set the logger to print only error messages</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">msLogger</span><span class="o">.</span><span class="n">setPriorityMin</span><span class="p">(</span><span class="n">Priority</span><span class="o">.</span><span class="n">ERROR</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot; Create a units manager object, it will be very often used to instantiate </span>
<span class="sd">the different classes: it defines the unit of the intput and output of the objects</span>
<span class="sd"> by default. A units manager object is also usefull for conversion purposes </span>
<span class="sd">(we&#39;ll see that later).</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">unit</span> <span class="o">=</span> <span class="n">msUnitsManager</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="s">&quot;Angstrom Degree amu kcal/mol ps&quot;</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot; Let&#39;s define the H2 ZMatrix:</span>
<span class="sd">H</span>
<span class="sd">H 0 0.74191</span>
<span class="sd">The z-matrix is the object who define the link between coordinates value</span>
<span class="sd">and position of the elements. If you want to define a custom problem </span>
<span class="sd">(e.g. a chain of pendulums), you only have to write a similar class and</span>
<span class="sd"> you should look at the msZMat source files for documentation.  </span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">zmat</span> <span class="o">=</span> <span class="n">msZMat</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="n">unit</span><span class="p">)</span><span class="o">.</span><span class="n">setId</span><span class="p">(</span><span class="s">&quot;zmat&quot;</span><span class="p">)</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="s">&quot;H H 0 0.74191&quot;</span><span class="p">)</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="sd">When you print an object, the all class hierarchy of the object is described.</span>
<span class="sd">At the first level (i.e. base class), you have the description of the &#39;database mapper&#39;:</span>
<span class="sd">its purpose is to translate the c++ class structure into tree structure in order to</span>
<span class="sd">be saved or loaded. </span>
<span class="sd">Under it, the first derived class from the tree mapper, the description of the </span>
<span class="sd">&#39;Physical interface&#39; interface. The purpose of this interface is to define a unique </span>
<span class="sd">framework to define units and parameters.</span>
<span class="sd">At the third level you have the description of the entity class: a collection </span>
<span class="sd">of elements with their position and masse.</span>
<span class="sd">Finally the description of the Z-matrix is providen.</span>
<span class="sd">You can obtain more information about all classes in their doxygen documentation.</span>
<span class="sd">&quot;&quot;&quot;</span> 
<span class="k">print</span> <span class="n">zmat</span>

<span class="sd">&quot;&quot;&quot; Now a system is defined. A system is just a collection of entity (parent class</span>
<span class="sd">of Z-Matrix), that you add to the container:</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">system</span> <span class="o">=</span> <span class="n">msSystem</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="n">unit</span><span class="p">)</span><span class="o">.</span><span class="n">addEntity</span><span class="p">(</span><span class="n">zmat</span><span class="p">)</span>
<span class="k">print</span> <span class="n">system</span>

<span class="sd">&quot;&quot;&quot; If you have VTK and if you linked atomism with it, you can display the system</span>
<span class="sd">by uncommenting the two next lines</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="c"># viewer = msVisualizerAtoms.New()</span>
<span class="c"># viewer.watch(system)</span>

<span class="sd">&quot;&quot;&quot;</span>
<span class="sd">At this point the system is defined: there is some elements (atoms here) associated</span>
<span class="sd">to a mass and their position is defined from some coordinates (Z-Matrix here). </span>
<span class="sd">The next step is the definition of motion: how the elements are moving w/ some </span>
<span class="sd">generalized coordinates, and what is the energy cost of it.</span>
<span class="sd">Lets consider the vibration H-H in H2.</span>
<span class="sd">step 1: definition of the generalized coordinates, it contains one variable,</span>
<span class="sd"> the separation H-H.</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">gcoors</span> <span class="o">=</span> <span class="n">msGeneralizedCoordinates</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="n">unit</span><span class="p">)</span>
<span class="n">q0</span>     <span class="o">=</span> <span class="n">msScalarVariable</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="s">&quot;Angstrom&quot;</span><span class="p">)</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="mf">0.7414</span><span class="p">,</span> <span class="mf">0.1</span><span class="p">,</span> <span class="mi">5</span><span class="p">,</span> <span class="mf">0.1</span><span class="p">,</span> <span class="mf">0.1</span><span class="p">)</span><span class="o">.</span><span class="n">setId</span><span class="p">(</span><span class="s">&quot;q0&quot;</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot; try &#39;q0.set?&#39; for the description of the method. </span>
<span class="sd">The generalized coordinates object is a container for scalar variable: </span>
<span class="sd">&quot;&quot;&quot;</span> 
<span class="n">gcoors</span><span class="o">.</span><span class="n">addVar</span><span class="p">(</span><span class="n">q0</span><span class="p">)</span>
<span class="k">print</span> <span class="n">gcoors</span>

<span class="sd">&quot;&quot;&quot; step 2: definition of the potential energy.</span>
<span class="sd">The potential energy can be any object deriving from scalar functions,</span>
<span class="sd">here we are using a parser type to define the Morse potential used by Martson </span>
<span class="sd"> and Bali-Kunti in their study. </span>
<span class="sd">see C. Clay Martson and G. G. Balint-Kurti, J Chem Phys 91(6): 3571 (1989)</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">epot</span> <span class="o">=</span> <span class="n">msParser</span><span class="o">.</span><span class="n">New</span><span class="p">(</span> <span class="n">unit</span><span class="o">.</span><span class="n">getEnergyStr</span><span class="p">()</span> <span class="p">)</span><span class="o">.</span><span class="n">setCoordinates</span><span class="p">(</span><span class="n">gcoors</span><span class="p">)</span>
<span class="n">epot</span><span class="o">.</span><span class="n">setExpression</span><span class="p">(</span><span class="s">&quot;De * ( 1 - exp( - beta*( q0 - re ) ))^2&quot;</span><span class="p">)</span><span class="o">.</span><span class="n">setId</span><span class="p">(</span><span class="s">&quot;f=Morse&quot;</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot; The unit object allow easy conversion of the tabulated parameters to </span>
<span class="sd">our unit system:</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">De</span>   <span class="o">=</span> <span class="n">unit</span><span class="o">.</span><span class="n">convert</span><span class="p">(</span><span class="s">&quot;eV&quot;</span><span class="p">,</span><span class="mf">4.7457</span><span class="p">)</span>
<span class="n">beta</span> <span class="o">=</span> <span class="n">unit</span><span class="o">.</span><span class="n">convert</span><span class="p">(</span><span class="s">&quot;m^-1&quot;</span><span class="p">,</span><span class="mf">1.94196e10</span><span class="p">)</span>
<span class="n">re</span>   <span class="o">=</span> <span class="n">unit</span><span class="o">.</span><span class="n">convert</span><span class="p">(</span><span class="s">&quot;m&quot;</span><span class="p">,</span><span class="mf">0.74191e-10</span><span class="p">)</span>

<span class="n">epot</span><span class="o">.</span><span class="n">setConstant</span><span class="p">(</span><span class="s">&quot;De&quot;</span><span class="p">,</span><span class="n">De</span><span class="p">)</span><span class="o">.</span><span class="n">setConstant</span><span class="p">(</span><span class="s">&quot;beta&quot;</span><span class="p">,</span><span class="n">beta</span><span class="p">)</span><span class="o">.</span><span class="n">setConstant</span><span class="p">(</span><span class="s">&quot;re&quot;</span><span class="p">,</span><span class="n">re</span><span class="p">)</span>
<span class="k">print</span> <span class="n">epot</span>

<span class="sd">&quot;&quot;&quot; step 3: definition of the kinetic operator. </span>
<span class="sd">The kinetic operator describes how the masses in the &#39;system&#39; object move with respect </span>
<span class="sd">to the generalized coordinates.</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">kinop</span> <span class="o">=</span> <span class="n">msKineticOperator</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="n">unit</span><span class="p">)</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="n">gcoors</span><span class="p">,</span><span class="n">system</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot; The kinetic operator is a container of dynamic Degrees Of Freedom (DOF);</span>
<span class="sd">these DOF link a coordinate in an entity to a scalar function.</span>
<span class="sd">Here, we want to link the separation H-H in the z-matrix (DOF number 0 of &#39;zmat&#39;) </span>
<span class="sd">by the scalar function f(q0)=q0.</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">kinop</span><span class="o">.</span><span class="n">addDynamicDof</span><span class="p">(</span> <span class="n">zmat</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="n">msParser</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="s">&quot;Angstrom&quot;</span><span class="p">)</span><span class="o">.</span><span class="n">setCoordinates</span><span class="p">(</span><span class="n">gcoors</span><span class="p">)</span><span class="o">.</span><span class="n">setExpression</span><span class="p">(</span><span class="s">&quot;q0&quot;</span><span class="p">)</span><span class="o">.</span><span class="n">setId</span><span class="p">(</span><span class="s">&quot;f=q0&quot;</span><span class="p">)</span> <span class="p">)</span>
<span class="k">print</span> <span class="n">kinop</span>

<span class="sd">&quot;&quot;&quot; PES and kinetic operator defined a motion, in the case of 1D motion, a</span>
<span class="sd">quantum resolution of the equation of motion is available.</span>
<span class="sd">&quot;&quot;&quot;</span> 
<span class="n">qmotion</span> <span class="o">=</span> <span class="n">msQuantumMotion1d</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="n">unit</span><span class="p">)</span><span class="o">.</span><span class="n">setEqOfMotion</span><span class="p">(</span><span class="n">kinop</span><span class="p">,</span><span class="n">epot</span><span class="p">)</span>
<span class="k">print</span> <span class="n">qmotion</span>

<span class="sd">&quot;&quot;&quot; In the previous ouput, you can look at the default value of the parameters.</span>
<span class="sd">To change the value of a parameter: </span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">qmotion</span><span class="o">.</span><span class="n">getParameters</span><span class="p">()</span><span class="o">.</span><span class="n">setParameter</span><span class="p">(</span><span class="s">&quot;BasisSize&quot;</span><span class="p">,</span><span class="mi">501</span><span class="p">,</span><span class="n">unit</span><span class="p">)</span>
<span class="n">qmotion</span><span class="o">.</span><span class="n">getParameters</span><span class="p">()</span><span class="o">.</span><span class="n">setParameter</span><span class="p">(</span><span class="s">&quot;Emax&quot;</span><span class="p">,</span><span class="n">De</span><span class="p">,</span><span class="n">unit</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot; All objects derived from the msMotion class contains several themodynamic and DOS related</span>
<span class="sd">functions; see the documentation of the class msMotion. For instance:</span>
<span class="sd">&quot;&quot;&quot;</span> 
<span class="n">Q</span>   <span class="o">=</span> <span class="n">qmotion</span><span class="o">.</span><span class="n">Q</span><span class="p">(</span><span class="mi">1000</span><span class="p">)</span>
<span class="n">S</span>   <span class="o">=</span> <span class="n">qmotion</span><span class="o">.</span><span class="n">S</span><span class="p">(</span><span class="mi">1000</span><span class="p">)</span>
<span class="n">Cp</span>  <span class="o">=</span> <span class="n">qmotion</span><span class="o">.</span><span class="n">Cp</span><span class="p">(</span><span class="mi">1000</span><span class="p">)</span>
<span class="n">dos</span> <span class="o">=</span> <span class="n">qmotion</span><span class="o">.</span><span class="n">DOS</span><span class="p">(</span><span class="mi">10</span><span class="p">)</span>
<span class="n">sos</span> <span class="o">=</span> <span class="n">qmotion</span><span class="o">.</span><span class="n">SOS</span><span class="p">(</span><span class="mi">10</span><span class="p">)</span>
<span class="sd">&quot;&quot;&quot; We want to compare the results with those of an harmonic oscilator.</span>
<span class="sd">From the PES &#39;epot&#39;, we define its harmonic approximation at the </span>
<span class="sd">equilibrium separation.</span>
<span class="sd">&quot;&quot;&quot;</span> 
<span class="n">epotHO</span> <span class="o">=</span> <span class="n">msParser</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="n">unit</span><span class="o">.</span><span class="n">getEnergyStr</span><span class="p">())</span><span class="o">.</span><span class="n">setCoordinates</span><span class="p">(</span><span class="n">gcoors</span><span class="p">)</span><span class="o">.</span><span class="n">setExpression</span><span class="p">(</span><span class="s">&quot;alpha * ( q0 - re )^2&quot;</span><span class="p">)</span><span class="o">.</span><span class="n">setId</span><span class="p">(</span><span class="s">&quot;f=HO&quot;</span><span class="p">)</span>

<span class="n">gcoors</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">Value</span> <span class="o">=</span> <span class="n">re</span>

<span class="sd">&quot;&quot;&quot; By default scalar functions implement a numeric derivation procedure.</span>
<span class="sd">The indice given in the function call is the indice of the variable we are deriving from.</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">epotHO</span><span class="o">.</span><span class="n">setConstant</span><span class="p">(</span><span class="s">&quot;alpha&quot;</span><span class="p">,</span> <span class="n">epot</span><span class="o">.</span><span class="n">partialSecondDerivative</span><span class="p">(</span><span class="mi">0</span><span class="p">)</span> <span class="p">)</span><span class="o">.</span><span class="n">setConstant</span><span class="p">(</span><span class="s">&quot;re&quot;</span><span class="p">,</span><span class="n">re</span><span class="p">)</span>

<span class="n">oscillsHO</span> <span class="o">=</span> <span class="n">msQuantumMotion1d</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="n">unit</span><span class="p">)</span><span class="o">.</span><span class="n">setEqOfMotion</span><span class="p">(</span><span class="n">kinop</span><span class="p">,</span><span class="n">epotHO</span><span class="p">)</span>
<span class="n">oscillsHO</span><span class="o">.</span><span class="n">getParameters</span><span class="p">()</span><span class="o">.</span><span class="n">setParameter</span><span class="p">(</span><span class="s">&quot;BasisSize&quot;</span><span class="p">,</span><span class="mi">501</span><span class="p">,</span><span class="s">&quot;&quot;</span><span class="p">)</span>
<span class="n">oscillsHO</span><span class="o">.</span><span class="n">getParameters</span><span class="p">()</span><span class="o">.</span><span class="n">setParameter</span><span class="p">(</span><span class="s">&quot;Emax&quot;</span><span class="p">,</span><span class="n">De</span><span class="p">,</span><span class="n">unit</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot; Now we need to define overall rotational and overall translational motions.</span>
<span class="sd">These modes of motion are not studied bt their equation of motion, but by </span>
<span class="sd">analytically derived statistical expressions.</span>
<span class="sd">Be carefull at the geometry of your system at the time your are fixing the </span>
<span class="sd">moment of inertia. Here, we recompute the atomic position at the current Value </span>
<span class="sd">of q0 (i.e. re, see 4 lines ahead)</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">kinop</span><span class="o">.</span><span class="n">setDynamicDofs</span><span class="p">()</span>
<span class="n">system</span><span class="o">.</span><span class="n">computeCartCoordinates</span><span class="p">()</span>
<span class="n">rotor</span> <span class="o">=</span> <span class="n">msLinearRigidRotor</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="n">unit</span><span class="p">)</span><span class="o">.</span><span class="n">setMmtOfInertia</span><span class="p">(</span><span class="n">system</span><span class="p">)</span><span class="o">.</span><span class="n">setSymmetryNumber</span><span class="p">(</span><span class="mi">2</span><span class="p">)</span>
<span class="n">rotor</span><span class="o">.</span><span class="n">getParameters</span><span class="p">()</span><span class="o">.</span><span class="n">setParameter</span><span class="p">(</span><span class="s">&quot;Emax&quot;</span><span class="p">,</span><span class="n">unit</span><span class="o">.</span><span class="n">convert</span><span class="p">(</span><span class="s">&quot;kcal/mol&quot;</span><span class="p">,</span><span class="mi">200</span><span class="p">),</span><span class="n">unit</span><span class="p">)</span>

<span class="n">trans</span> <span class="o">=</span> <span class="n">msTranslation</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="n">unit</span><span class="p">)</span><span class="o">.</span><span class="n">setSystem</span><span class="p">(</span><span class="n">system</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot; At this point everything is defined. Let&#39;s realized some computation</span>
<span class="sd">and plot the results.</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">qmotion</span><span class="o">.</span><span class="n">computeEigenValues</span><span class="p">()</span>
<span class="n">oscillsHO</span><span class="o">.</span><span class="n">computeEigenValues</span><span class="p">()</span>

<span class="sd">&quot;&quot;&quot; In this particular case the kinetic function (or reduced mass because q0 is a length)</span>
<span class="sd"> is constant w/q0. That means that the kinetic function can be freezed to its current </span>
<span class="sd">value.</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">qmotion</span><span class="o">.</span><span class="n">freezeKinValue</span><span class="p">(</span> <span class="n">kinop</span><span class="o">.</span><span class="n">kineticFunction</span><span class="p">()</span> <span class="p">)</span>
<span class="n">oscillsHO</span><span class="o">.</span><span class="n">freezeKinValue</span><span class="p">(</span> <span class="n">kinop</span><span class="o">.</span><span class="n">kineticFunction</span><span class="p">()</span> <span class="p">)</span>

<span class="sd">&quot;&quot;&quot; get the number 0, 4 and 9 eigenvectors, they will be plot later </span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">ev0</span>  <span class="o">=</span> <span class="n">qmotion</span><span class="o">.</span><span class="n">getEigenvector</span><span class="p">(</span><span class="mi">0</span><span class="p">)</span>
<span class="n">ev5</span>  <span class="o">=</span> <span class="n">qmotion</span><span class="o">.</span><span class="n">getEigenvector</span><span class="p">(</span><span class="mi">4</span><span class="p">)</span>
<span class="n">ev10</span> <span class="o">=</span> <span class="n">qmotion</span><span class="o">.</span><span class="n">getEigenvector</span><span class="p">(</span><span class="mi">9</span><span class="p">)</span>
<span class="n">ev0HO</span> <span class="o">=</span> <span class="n">oscillsHO</span><span class="o">.</span><span class="n">getEigenvector</span><span class="p">(</span><span class="mi">0</span><span class="p">)</span>
<span class="n">ev5HO</span>  <span class="o">=</span> <span class="n">oscillsHO</span><span class="o">.</span><span class="n">getEigenvector</span><span class="p">(</span><span class="mi">4</span><span class="p">)</span>
<span class="n">ev10HO</span> <span class="o">=</span> <span class="n">oscillsHO</span><span class="o">.</span><span class="n">getEigenvector</span><span class="p">(</span><span class="mi">9</span><span class="p">)</span>


<span class="n">rotor</span><span class="o">.</span><span class="n">computeEigenValues</span><span class="p">()</span>
<span class="sd">&quot;&quot;&quot; Occupency probability of the rotor eigenmodes at 500 and 2000 K</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">Jdist500</span><span class="o">=</span><span class="n">rotor</span><span class="o">.</span><span class="n">getOccupencyProbability</span><span class="p">(</span><span class="mi">500</span><span class="p">)</span>
<span class="n">Jdist2000</span><span class="o">=</span><span class="n">rotor</span><span class="o">.</span><span class="n">getOccupencyProbability</span><span class="p">(</span><span class="mi">2000</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot; Compute the total entropy for the Morse and Harmonic PES </span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">s</span><span class="o">=</span><span class="p">[]</span>
<span class="n">sHO</span><span class="o">=</span><span class="p">[]</span>
<span class="n">Tlist</span>  <span class="o">=</span> <span class="p">[</span><span class="mi">200</span><span class="p">,</span> <span class="mi">300</span><span class="p">,</span> <span class="mi">500</span><span class="p">,</span> <span class="mi">800</span><span class="p">,</span> <span class="mi">1200</span><span class="p">,</span> <span class="mi">1500</span><span class="p">,</span> <span class="mi">2000</span><span class="p">,</span> <span class="mi">2500</span><span class="p">,</span> <span class="mi">3000</span><span class="p">]</span> 
<span class="k">for</span> <span class="n">T</span> <span class="ow">in</span> <span class="n">Tlist</span><span class="p">:</span> 
    <span class="n">s</span><span class="o">.</span><span class="n">append</span><span class="p">(</span>   <span class="n">trans</span><span class="o">.</span><span class="n">S</span><span class="p">(</span><span class="n">T</span><span class="p">)</span>  <span class="o">+</span> <span class="n">rotor</span><span class="o">.</span><span class="n">S</span><span class="p">(</span><span class="n">T</span><span class="p">)</span> <span class="o">+</span> <span class="n">qmotion</span><span class="o">.</span><span class="n">S</span><span class="p">(</span><span class="n">T</span><span class="p">)</span>   <span class="o">+</span> <span class="n">unit</span><span class="o">.</span><span class="n">convert</span><span class="p">(</span><span class="s">&quot;J/mol&quot;</span><span class="p">,</span><span class="n">R</span><span class="p">)</span> <span class="p">)</span>
    <span class="n">sHO</span><span class="o">.</span><span class="n">append</span><span class="p">(</span> <span class="n">trans</span><span class="o">.</span><span class="n">S</span><span class="p">(</span><span class="n">T</span><span class="p">)</span>  <span class="o">+</span> <span class="n">rotor</span><span class="o">.</span><span class="n">S</span><span class="p">(</span><span class="n">T</span><span class="p">)</span> <span class="o">+</span> <span class="n">oscillsHO</span><span class="o">.</span><span class="n">S</span><span class="p">(</span><span class="n">T</span><span class="p">)</span> <span class="o">+</span> <span class="n">unit</span><span class="o">.</span><span class="n">convert</span><span class="p">(</span><span class="s">&quot;J/mol&quot;</span><span class="p">,</span><span class="n">R</span><span class="p">)</span> <span class="p">)</span>

<span class="sd">&quot;&quot;&quot; Same for calorific capacity</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">cp</span><span class="o">=</span><span class="p">[]</span>
<span class="n">cpHO</span><span class="o">=</span><span class="p">[]</span>
<span class="k">for</span> <span class="n">T</span> <span class="ow">in</span> <span class="n">Tlist</span><span class="p">:</span>
    <span class="n">cp</span><span class="o">.</span><span class="n">append</span><span class="p">(</span>   <span class="n">trans</span><span class="o">.</span><span class="n">Cv</span><span class="p">(</span><span class="n">T</span><span class="p">)</span> <span class="o">+</span> <span class="n">rotor</span><span class="o">.</span><span class="n">Cv</span><span class="p">(</span><span class="n">T</span><span class="p">)</span> <span class="o">+</span> <span class="n">qmotion</span><span class="o">.</span><span class="n">Cv</span><span class="p">(</span><span class="n">T</span><span class="p">)</span>   <span class="o">+</span>  <span class="n">unit</span><span class="o">.</span><span class="n">convert</span><span class="p">(</span><span class="s">&quot;J/mol&quot;</span><span class="p">,</span><span class="n">R</span><span class="p">)</span> <span class="p">)</span>
    <span class="n">cpHO</span><span class="o">.</span><span class="n">append</span><span class="p">(</span> <span class="n">trans</span><span class="o">.</span><span class="n">Cv</span><span class="p">(</span><span class="n">T</span><span class="p">)</span> <span class="o">+</span> <span class="n">rotor</span><span class="o">.</span><span class="n">Cv</span><span class="p">(</span><span class="n">T</span><span class="p">)</span> <span class="o">+</span> <span class="n">oscillsHO</span><span class="o">.</span><span class="n">Cv</span><span class="p">(</span><span class="n">T</span><span class="p">)</span> <span class="o">+</span>  <span class="n">unit</span><span class="o">.</span><span class="n">convert</span><span class="p">(</span><span class="s">&quot;J/mol&quot;</span><span class="p">,</span><span class="n">R</span><span class="p">)</span> <span class="p">)</span>

    
<span class="sd">&quot;&quot;&quot; See the documentation of matplotlib for the plotting part.</span>
<span class="sd">All the object derived from a scalar function can return a python list to </span>
<span class="sd">be plotted using the &#39;getXSample1d&#39; and &#39;getYSample1d&#39; methods. </span>
<span class="sd">In these functions the first variable ie the indice of the coordinates </span>
<span class="sd">along witch you want to sample, and the second the number of sample you want.</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">fig</span><span class="p">,</span> <span class="n">axes</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">subplots</span><span class="p">(</span><span class="mi">3</span><span class="p">,</span> <span class="mi">2</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot; The following four lines define the analytical eigenvalues of the Morse potential </span>
<span class="sd">in the H2 molecule. see C. Clay Martson and G. G. Balint-Kurti, J Chem Phys 91(6): 3571 (1989)</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">scoor</span><span class="o">=</span><span class="p">[</span><span class="n">q0</span><span class="o">.</span><span class="n">Min</span><span class="p">,</span><span class="n">q0</span><span class="o">.</span><span class="n">Max</span><span class="p">]</span>
<span class="n">y0_analytic</span><span class="o">=</span><span class="p">[</span><span class="n">unit</span><span class="o">.</span><span class="n">convert</span><span class="p">(</span><span class="s">&quot;hartree&quot;</span><span class="p">,</span><span class="mf">9.86922e-3</span><span class="p">),</span><span class="n">unit</span><span class="o">.</span><span class="n">convert</span><span class="p">(</span><span class="s">&quot;hartree&quot;</span><span class="p">,</span>  <span class="mf">9.86922e-3</span><span class="p">)]</span>
<span class="n">y4_analytic</span><span class="o">=</span><span class="p">[</span><span class="n">unit</span><span class="o">.</span><span class="n">convert</span><span class="p">(</span><span class="s">&quot;hartree&quot;</span><span class="p">,</span><span class="mf">7.847518e-2</span><span class="p">),</span><span class="n">unit</span><span class="o">.</span><span class="n">convert</span><span class="p">(</span><span class="s">&quot;hartree&quot;</span><span class="p">,</span> <span class="mf">7.847518e-2</span><span class="p">)]</span>
<span class="n">y9_analytic</span><span class="o">=</span><span class="p">[</span><span class="n">unit</span><span class="o">.</span><span class="n">convert</span><span class="p">(</span><span class="s">&quot;hartree&quot;</span><span class="p">,</span><span class="mf">1.3836303e-1</span><span class="p">),</span><span class="n">unit</span><span class="o">.</span><span class="n">convert</span><span class="p">(</span><span class="s">&quot;hartree&quot;</span><span class="p">,</span><span class="mf">1.3836303e-1</span><span class="p">)]</span>

<span class="sd">&quot;&quot;&quot;on top right subplot, we plot the potential energy and the eigenvectors with </span>
<span class="sd"> a sampling of 500 along the first coordinates of &#39;gencoors&#39; (indiced 0, i.e. q0);</span>
<span class="sd"> the number &#39;600&#39; is a scaling factor for having a nice drawing of the eigenvectors</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span>    <span class="n">epotHO</span><span class="o">.</span><span class="n">getXSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">),</span> <span class="n">epotHO</span><span class="o">.</span><span class="n">getYSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">),</span> <span class="s">&quot;gray&quot;</span><span class="p">,</span> 
	            <span class="n">ev0HO</span><span class="o">.</span><span class="n">getXSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">)</span> <span class="p">,</span> <span class="p">[</span> <span class="mi">600</span><span class="o">*</span><span class="n">x</span> <span class="o">+</span> <span class="n">oscillsHO</span><span class="o">.</span><span class="n">eigenvalue</span><span class="p">(</span><span class="mi">0</span><span class="p">)</span>  <span class="k">for</span> <span class="n">x</span> <span class="ow">in</span> <span class="n">ev0HO</span><span class="o">.</span><span class="n">getYSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">)]</span> <span class="p">,</span><span class="s">&quot;gray&quot;</span><span class="p">,</span>  
	            <span class="n">ev5HO</span><span class="o">.</span><span class="n">getXSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">)</span> <span class="p">,</span> <span class="p">[</span> <span class="mi">600</span><span class="o">*</span><span class="n">x</span> <span class="o">+</span> <span class="n">oscillsHO</span><span class="o">.</span><span class="n">eigenvalue</span><span class="p">(</span><span class="mi">4</span><span class="p">)</span>  <span class="k">for</span> <span class="n">x</span> <span class="ow">in</span> <span class="n">ev5HO</span><span class="o">.</span><span class="n">getYSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">)]</span> <span class="p">,</span><span class="s">&quot;gray&quot;</span><span class="p">,</span>
	            <span class="n">ev10HO</span><span class="o">.</span><span class="n">getXSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">),</span> <span class="p">[</span> <span class="mi">600</span><span class="o">*</span><span class="n">x</span> <span class="o">+</span> <span class="n">oscillsHO</span><span class="o">.</span><span class="n">eigenvalue</span><span class="p">(</span><span class="mi">9</span><span class="p">)</span> <span class="k">for</span> <span class="n">x</span> <span class="ow">in</span> <span class="n">ev10HO</span><span class="o">.</span><span class="n">getYSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">)],</span><span class="s">&quot;gray&quot;</span><span class="p">,</span>
		    <span class="n">epot</span><span class="o">.</span><span class="n">getXSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">),</span> <span class="n">epot</span><span class="o">.</span><span class="n">getYSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">),</span> <span class="s">&quot;blue&quot;</span><span class="p">,</span> 
	            <span class="n">ev0</span><span class="o">.</span><span class="n">getXSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">)</span> <span class="p">,</span> <span class="p">[</span> <span class="mi">600</span><span class="o">*</span><span class="n">x</span> <span class="o">+</span> <span class="n">qmotion</span><span class="o">.</span><span class="n">eigenvalue</span><span class="p">(</span><span class="mi">0</span><span class="p">)</span>  <span class="k">for</span> <span class="n">x</span> <span class="ow">in</span> <span class="n">ev0</span><span class="o">.</span><span class="n">getYSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">)]</span> <span class="p">,</span><span class="s">&quot;blue&quot;</span><span class="p">,</span>  
	            <span class="n">ev5</span><span class="o">.</span><span class="n">getXSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">)</span> <span class="p">,</span> <span class="p">[</span> <span class="mi">600</span><span class="o">*</span><span class="n">x</span> <span class="o">+</span> <span class="n">qmotion</span><span class="o">.</span><span class="n">eigenvalue</span><span class="p">(</span><span class="mi">4</span><span class="p">)</span>  <span class="k">for</span> <span class="n">x</span> <span class="ow">in</span> <span class="n">ev5</span><span class="o">.</span><span class="n">getYSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">)]</span> <span class="p">,</span><span class="s">&quot;blue&quot;</span><span class="p">,</span>
	            <span class="n">ev10</span><span class="o">.</span><span class="n">getXSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">),</span> <span class="p">[</span> <span class="mi">600</span><span class="o">*</span><span class="n">x</span> <span class="o">+</span> <span class="n">qmotion</span><span class="o">.</span><span class="n">eigenvalue</span><span class="p">(</span><span class="mi">9</span><span class="p">)</span> <span class="k">for</span> <span class="n">x</span> <span class="ow">in</span> <span class="n">ev10</span><span class="o">.</span><span class="n">getYSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">)],</span><span class="s">&quot;blue&quot;</span><span class="p">,</span>
		    <span class="n">scoor</span><span class="p">,</span><span class="n">y0_analytic</span><span class="p">,</span><span class="s">&quot;red&quot;</span><span class="p">,</span>
		    <span class="n">scoor</span><span class="p">,</span><span class="n">y4_analytic</span><span class="p">,</span><span class="s">&quot;red&quot;</span><span class="p">,</span>
	            <span class="n">scoor</span><span class="p">,</span><span class="n">y9_analytic</span><span class="p">,</span><span class="s">&quot;red&quot;</span><span class="p">,</span>
		    <span class="p">)</span>
	            
<span class="n">axes</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">axis</span><span class="p">((</span><span class="mi">0</span><span class="p">,</span><span class="mf">3.5</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="n">unit</span><span class="o">.</span><span class="n">convert</span><span class="p">(</span><span class="s">&quot;kcal/mol&quot;</span><span class="p">,</span><span class="mi">160</span><span class="p">)))</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s">&quot;H-H [&quot;</span><span class="o">+</span><span class="n">unit</span><span class="o">.</span><span class="n">getLengthStr</span><span class="p">()</span><span class="o">+</span><span class="s">&quot;]&quot;</span><span class="p">,</span><span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s">&quot;energy [&quot;</span><span class="o">+</span><span class="n">unit</span><span class="o">.</span><span class="n">getEnergyStr</span><span class="p">()</span><span class="o">+</span><span class="s">&quot;]&quot;</span><span class="p">,</span><span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">tick_params</span><span class="p">(</span><span class="n">labelsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>

<span class="c"># On the top left subplot we plot the sum of states of &#39;qmotion&#39; and compare it </span>
<span class="c"># with the one of the harmonic oscillator:</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span> <span class="n">qmotion</span><span class="o">.</span><span class="n">SOS</span><span class="p">()</span><span class="o">.</span><span class="n">getXSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">),</span>    <span class="n">qmotion</span><span class="o">.</span><span class="n">SOS</span><span class="p">()</span><span class="o">.</span><span class="n">getYSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">),</span> <span class="s">&quot;-&quot;</span><span class="p">,</span>
                 <span class="n">oscillsHO</span><span class="o">.</span><span class="n">SOS</span><span class="p">()</span><span class="o">.</span><span class="n">getXSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">),</span>  <span class="n">oscillsHO</span><span class="o">.</span><span class="n">SOS</span><span class="p">()</span><span class="o">.</span><span class="n">getYSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">500</span><span class="p">),</span> <span class="s">&quot;gray&quot;</span><span class="p">)</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s">&quot;vibrational energy [&quot;</span><span class="o">+</span><span class="n">unit</span><span class="o">.</span><span class="n">getEnergyStr</span><span class="p">()</span><span class="o">+</span><span class="s">&quot;]&quot;</span><span class="p">,</span><span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s">&quot;vibrational sum of states&quot;</span><span class="p">,</span><span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">tick_params</span><span class="p">(</span><span class="n">labelsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
	
<span class="c"># On the middle left subplot we plot the occupency probabilities of the rotor w/</span>
<span class="c"># to the eigenvalues:</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span> <span class="n">Jdist500</span><span class="o">.</span><span class="n">getXSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">1000</span><span class="p">),</span>   <span class="n">Jdist500</span><span class="o">.</span><span class="n">getYSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">1000</span><span class="p">),</span>  <span class="s">&quot;-&quot;</span><span class="p">,</span>
                 <span class="n">Jdist2000</span><span class="o">.</span><span class="n">getXSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">1000</span><span class="p">),</span>  <span class="n">Jdist2000</span><span class="o">.</span><span class="n">getYSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">1000</span><span class="p">),</span> <span class="s">&quot;-&quot;</span><span class="p">)</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s">&quot;angular momentum number&quot;</span><span class="p">,</span><span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s">&quot;Occupency probability&quot;</span><span class="p">,</span><span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">tick_params</span><span class="p">(</span><span class="n">labelsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>

<span class="c"># On the middle right sublot we plot the rotor Density Of States and Sum Of States,</span>
<span class="c"># We need in a high resolution to represent the diracs of the DOS</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">rotor</span><span class="o">.</span><span class="n">DOS</span><span class="p">()</span><span class="o">.</span><span class="n">getXSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">10000</span><span class="p">),</span>  <span class="n">rotor</span><span class="o">.</span><span class="n">DOS</span><span class="p">()</span><span class="o">.</span><span class="n">getYSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">10000</span><span class="p">),</span> <span class="s">&quot;-&quot;</span><span class="p">,</span>
                <span class="n">rotor</span><span class="o">.</span><span class="n">SOS</span><span class="p">()</span><span class="o">.</span><span class="n">getXSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">1000</span><span class="p">),</span>   <span class="n">rotor</span><span class="o">.</span><span class="n">SOS</span><span class="p">()</span><span class="o">.</span><span class="n">getYSample1d</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">1000</span><span class="p">),</span> <span class="s">&quot;-&quot;</span><span class="p">)</span>  
<span class="n">axes</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s">&quot;rotational energy [&quot;</span><span class="o">+</span><span class="n">unit</span><span class="o">.</span><span class="n">getEnergyStr</span><span class="p">()</span><span class="o">+</span><span class="s">&quot;]&quot;</span><span class="p">,</span><span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s">&quot;rotational SOS and DOS&quot;</span><span class="p">,</span><span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">tick_params</span><span class="p">(</span><span class="n">labelsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">axis</span><span class="p">((</span><span class="mi">0</span><span class="p">,</span><span class="n">unit</span><span class="o">.</span><span class="n">convert</span><span class="p">(</span><span class="s">&quot;kcal/mol&quot;</span><span class="p">,</span><span class="mi">200</span><span class="p">),</span><span class="mi">0</span><span class="p">,</span><span class="mi">1400</span><span class="p">))</span>
<span class="c"># Now we compute the entropy and calorific capacity lists and plot them in a subplots.</span>
<span class="c"># Again, the results are compared to those of the harmonic oscillator.</span>

<span class="n">axes</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span> <span class="n">Tlist</span><span class="p">,</span> <span class="n">s</span> <span class="p">,</span> <span class="s">&quot;-&quot;</span><span class="p">,</span> <span class="n">Tlist</span><span class="p">,</span> <span class="n">sHO</span><span class="p">)</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s">&quot;Temperature [&quot;</span><span class="o">+</span><span class="n">unit</span><span class="o">.</span><span class="n">getTempStr</span><span class="p">()</span><span class="o">+</span><span class="s">&quot;]&quot;</span><span class="p">,</span><span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s">&quot;Total entropy [&quot;</span><span class="o">+</span><span class="n">unit</span><span class="o">.</span><span class="n">getEnergyStr</span><span class="p">()</span><span class="o">+</span><span class="s">&quot;]&quot;</span><span class="p">,</span><span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">tick_params</span><span class="p">(</span><span class="n">labelsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>

<span class="n">axes</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span> <span class="n">Tlist</span><span class="p">,</span> <span class="n">cp</span><span class="p">,</span> <span class="s">&quot;-&quot;</span><span class="p">,</span> <span class="n">Tlist</span><span class="p">,</span> <span class="n">cpHO</span><span class="p">,</span><span class="s">&#39;gray&#39;</span><span class="p">)</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s">&quot;Temperature [&quot;</span><span class="o">+</span><span class="n">unit</span><span class="o">.</span><span class="n">getTempStr</span><span class="p">()</span><span class="o">+</span><span class="s">&quot;]&quot;</span><span class="p">,</span><span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s">&quot;Total Cp [&quot;</span><span class="o">+</span><span class="n">unit</span><span class="o">.</span><span class="n">getEnergyStr</span><span class="p">()</span><span class="o">+</span><span class="s">&quot;]&quot;</span><span class="p">,</span><span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">tick_params</span><span class="p">(</span><span class="n">labelsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>

<span class="c"># this is the experimental values, coming from the janaf table,  </span>
<span class="c"># see http://kinetics.nist.gov/janaf/html/H-050.html, defined in J/mol/K</span>
<span class="n">TListexp</span>  <span class="o">=</span>  <span class="p">[</span><span class="mi">200</span><span class="p">,</span>     <span class="mi">300</span><span class="p">,</span>     <span class="mi">500</span><span class="p">,</span>     <span class="mi">800</span><span class="p">,</span>     <span class="mi">1200</span><span class="p">,</span>   <span class="mi">1500</span><span class="p">,</span>    <span class="mi">2000</span><span class="p">,</span>    <span class="mi">2500</span><span class="p">,</span>    <span class="mi">3000</span>   <span class="p">]</span> 
<span class="n">ExpListS</span>   <span class="o">=</span> <span class="p">[</span><span class="mf">119.412</span><span class="p">,</span> <span class="mf">130.858</span><span class="p">,</span> <span class="mf">145.737</span><span class="p">,</span> <span class="mf">159.548</span><span class="p">,</span> <span class="mf">171.79</span><span class="p">,</span> <span class="mf">178.846</span><span class="p">,</span> <span class="mf">188.418</span><span class="p">,</span> <span class="mf">196.243</span><span class="p">,</span> <span class="mf">202.891</span><span class="p">]</span>   
<span class="n">ExpListCp</span> <span class="o">=</span>  <span class="p">[</span><span class="mf">27.447</span><span class="p">,</span>  <span class="mf">28.849</span><span class="p">,</span>  <span class="mf">29.260</span><span class="p">,</span>  <span class="mf">29.624</span><span class="p">,</span>  <span class="mf">30.992</span><span class="p">,</span> <span class="mf">32.298</span><span class="p">,</span>  <span class="mf">34.280</span><span class="p">,</span>  <span class="mf">35.842</span><span class="p">,</span>  <span class="mf">37.087</span> <span class="p">]</span> 
 
<span class="n">axes</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span> <span class="n">TListexp</span><span class="p">,</span> <span class="p">[</span> <span class="n">unit</span><span class="o">.</span><span class="n">convert</span><span class="p">(</span><span class="s">&quot;J/mol&quot;</span><span class="p">,</span><span class="n">x</span><span class="p">)</span> <span class="k">for</span> <span class="n">x</span> <span class="ow">in</span> <span class="n">ExpListS</span><span class="p">],</span>   <span class="s">&quot;r*&quot;</span> <span class="p">,</span><span class="n">markersize</span><span class="o">=</span><span class="mi">15</span><span class="p">)</span>
<span class="n">axes</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span> <span class="n">TListexp</span><span class="p">,</span> <span class="p">[</span> <span class="n">unit</span><span class="o">.</span><span class="n">convert</span><span class="p">(</span><span class="s">&quot;J/mol&quot;</span><span class="p">,</span><span class="n">x</span><span class="p">)</span> <span class="k">for</span> <span class="n">x</span> <span class="ow">in</span> <span class="n">ExpListCp</span><span class="p">],</span>  <span class="s">&quot;r*&quot;</span> <span class="p">,</span><span class="n">markersize</span><span class="o">=</span><span class="mi">15</span><span class="p">)</span>

<span class="n">fig</span><span class="o">.</span><span class="n">set_size_inches</span><span class="p">(</span><span class="mi">20</span><span class="p">,</span><span class="mi">20</span><span class="p">)</span>
<span class="n">fig</span><span class="o">.</span><span class="n">savefig</span><span class="p">(</span><span class="s">&#39;H2_0.png&#39;</span><span class="p">)</span>
<span class="n">fig</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>


 
</pre></div>
</body>
</html>
